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This paper describes model structures and parameter estimation algorithms suit- 
able for the identification of unsteady aerodynamic models from input-output data. 
The model structures presented are state space models and include linear time- 
invariant (LTI) models and linear parameter-varying (LPV) models. They cover 
a wide range of local and parameter dependent identification problems arising in 
unsteady aerodynamics and nonlinear flight dynamics. We present a residue algo- 
rithm for estimating model parameters from data. The algorithm can incorporate 
apriori information and is described in detail. The algorithms are evaluated on the 
F-16XL wind-tunnel test data from NASA Langley Research Center. Results of 
numerical evaluation are presented. The paper concludes with a discussion major 
issues and directions for future work. 


I. Introduction 


A model is a useful representation of system dynamics. The precise meaning of “useful” is given 
by the application for which the model is intended. As a general rule, useful models of engineering 
systems must be computationally efficient. Another important requirement is the ability to ac- 
curately describe system behaviors under different conditions, i.e. , predictive power of the model. 
These two requirements are competing interests in that models that predict many system behaviors 
tend to be computationally hard and vice versa. For instance, Navier-Stokes equations have a lot 
more predictive power than the aerodynamic models discussed in this paper, but they are also 
numerically much harder. The modeler’s task is then to develop models that balance computa- 
tional and predictive requirements perhaps by neglecting fine scale structure and considering only 
those features that are dominant in the application. This is, of course, easier said than done in 
most flight dynamics problems where distributed phenomena interact to produce complex dynam- 
ical behaviors. A good example is tail-buffeting at high angles of attack in which vortices shed 
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from wing impinge tail surface causing structural vibrations and eventually fatigue. Identification 
of models from experimental data may be the only avenue in such problems. In other problems 
where first principles can be used to derive models such as in rigid body dynamics, identification 
of model parameters such as stability derivatives from experimental data is required. Researchers 
have also developed simple models for complicated but special flow problems by considering macro- 
scale flow physics. Such models also contain parameters that must be estimated from experimental 
data. 5,6,8,10 



Not Validated 


Figure 1. Steps in system identification 


System identification is an iterative process consisting of experiment design, model structure de- 
termination, parameter estimation and model validation (Figure l). 7,13 Experiment design is con- 
cerned with the selection of test inputs and measurement variables, their sampling and condition- 
ing, and other aspects that are directly related to the notion of an informative experiment. An 
experiment is said to be informative with respect to a model set if the experimental data allows 
discrimination between models in the set. Model structure selection and parameter estimation are 
simultaneous tasks, though in practice, they are performed sequentially with model structure being 
selected first. The structure of a model is defined in terms of its nature (eg. linear black-box, 
Wiener, neural-network), order (number of free variables) and parameterization (e.g state space in 
modal coordinates). There are several considerations that go into model structure selection includ- 
ing apriori knowledge, flexibility and algorithmic complexity. After selecting a model structure, 
parameters of the model are estimated from experimental data. Minimum mean square estima- 
tion (MMSE) and maximum likelihood estimation (MLE) 7, 13,15 are the most common methods 
of parameter estimation. Numerical procedures for estimation range from least squares to the 
expectation-maximization (EM) algorithm. The final step of model validation is to determine the 
“goodness” of an identified model, i.e, the extent to which it is useful in its intended application. 
Thus, a model identified for controller design may be deemed valid if the resulting controller pro- 
duces satisfactory performance. When model validation fails, it may be due to the choice of model 
structure and the lack of information in experimental data. As shown in Figure 1, the steps are 
iterated until a validated model is found. Details of these steps and the iterative process can be 
found in Ref. 13, 14, 15. 

This paper a describes time domain identification methods suitable for the modeling of unsteady 
and nonlinear aerodynamics behaviors arising in flight dynamics. Section II begins with model 

a A detailed version of the paper can be obtained from the authors 
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structures for modeling dynamics in the neighborhood of an operating condition, for example, fixed 
Mach number or fixed dynamic pressure. Such models are local models with a limited region of 
validity. The section also presents parameter-dependent model structures for modeling dynamics 
over a range of parameter values. We refer to such models as global models though they are not 
suitable models for certain global behaviors. The model structures are illustrated with examples 
relating them to well-known models such as the indicial response function models. A residue 
based algorithm for parameter estimation is given in Section III. The algorithms are evaluated 
on the F-16XL wind-tunnel test data from NASA Langley Research Center. Results of numerical 
evaluation are presented in Section IV. Conclusions and recommendations for future work are given 
in Section V. 

System identification methods have been used widely in aircraft applications. The special issue of 
the AIAA Journal of Aircraft 1 and the invited session at the AIAA Atmospheric Flight Mechanics 2 
contain papers describing many techniques and a broad range of recent applications. The main 
contributions of this paper are a new algorithm for identifying state space models of unsteady 
aerodynamic effects and its limited validation using F-16XL wind-tunnel test data. 

II. Local and Global Parameter Dependent Model Structures 


A model structure specifies how test inputs and external noises affect observed outputs. We only 
describe the relationship between test inputs and observed outputs to keep the presentation simple. 


A. Unstructured Linear Time-Invariant (LTI) Models 

Consider a continuous-time LTI state space model: 

x = Ax + Bu (la) 

y = Cx + Du (lb) 

where A € IR nX ”, B € IR nXm , C 6 Ai lyn , and D e IR Zxm . This model is completely determined 
by an order parameter n and the state space matrices (A. B. C. D). It is an unstructured model 
because the state space matrices have no specific structure and are allowed to take on any value. 
When the order parameter n is fixed, we get a fixed order unstructured model. 


B. Structured LTI Models 


In many applications, the natural frequencies and damping coefficients of system modes are known 
apriori with some confidence, e.g, cu n ±0.2u) n for a 20% variation from the nominal natural frequency 
co n . Physical reasoning may suggest the use of first order exponentially decaying or the second order 
decaying oscillatory model ,,6,8,1( or a combination of both. Thus, the A matrix takes the form: 


A = diagonal jr'i, 


Uv r , 


' 0 1 ' 
,-wf -2CiWi 


' 0 1 

.-^Nc —2 (n c Wn c _ 


( 2 ) 
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where {rk}j(fi are the real eigenvalues and {(u;*;, (k)}k=i are the eigenvalues with non-zero imaginary 
parts written in terms of natural frequency u & and damping coefficient 

The model given by (1) and (2) is specified by the order parameters ( N r ,N c ), the eigenvalue 
locations {rkf^Li and {(uik, Cfc)}£iJ and the matrices (B,C,D). If the postulated model (1) and 
(2) is used for system identification, some or all parameters in A may be known apriori with a 
given level of confidence. Then, apriori information on the parameters is described by a set A of 
admissible eigenvalue locations: 

A = {(<**, &)}£i) : rf° < Tj < r““, df* < < o£“, ■ ■ •} (3) 

where the bounds specify confidence intervals for each eigenvalue. 


Example II. 1 (Structured model) Consider an indicial response function model: 

y(t) — ci [ e ai ^ T ^a(T)dr + C 2 [ e a2 ^ T ^a(T)dT (4) 

Jo Jo 

= fGit-r) 

Jo 

where a\ < 0, 02 < 0, c\ and C 2 are unknown real numbers. This dynamical system in more general 
form appears in the work of Klein and Murphy 10 where it models the unsteady contribution to aero- 
dynamic coefficients in one degree of freedom motion. The strict negativity of a\ and a 2 guarantee 
that bounded test inputs (a and a) produce bounded outputs (y). The integral representation can 
be converted to state space form as: 



x = 


"ai 0 

.0 a- 2 _ 

y = [c 1 c- 2 ]x 


x + 


'1 O' 
0 1 


a 

a 


(6a) 

(6b) 


A relationship between model equations (6) and aircraft equations of motion will be shown in Sec- 
tion IV. In (6a), the A-matrix has a specific structure. Let us assume that a\ satisfies: 


Omin Ol f; Omax 0 


and that no information is available on the remaining parameters (other than strict negativity of 
a 2 ). Then, A takes the form: 


A — {( 01 , 02 ) : O m j n < Gi < O max <0, 02 < 0} 

Thus, indicial response function models of the form (4) are structured LTI models. We make two 
important observations. First, the number of (linearly independent) terms in the convolution form 
is equal to the model order parameter (state space dimension) . Second, the B-matrix is a constant. 
This is true of all LTI models. 
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C. Linear Parameter Varying (LPV) Global Models 


Many aerodynamic and flight dynamics behaviors result from changes in system dynamics brought 
about by changes in physical parameters. The LTI models described in the previous section are 
not appropriate for such systems. We present generalizations of the local models of the previous 
sections that can capture parameter-dependent phenomena in this section. 

Consider the system: 


x = A(a)x + B(a)u (7a) 

y = C(a)x + D(a)u (7b) 

where A(a) E TR nXn , B(a) E TR nXm , C(a) 6 IR 1XTi , and D(a) 6 IR 1Xm . A system of this form 
is called a linear parameter-varying (LPV) system. The system matrices are functions of angle of 
attack a and, during system operation, the angle of attack changes with time. Thus, along each 
possible a trajectory, (7) is a linear time-varying system. Since the specific trajectory taken by a 
can change from one run to another, the above system is not a single linear time-varying system, 
but a collection of linear time- varying systems. It is a global model in this restricted sense. 


Example II. 2 (Indicial function vs. LPV model) The integral kernel in Example II. 1 has the 

general form: 

G{t - t; a(r)) = J2G k {t- T)a(r) k 
k 

which implies that the dynamical part of those models, given by G k ’s, is LTI. Consider the state 
space model: 


x = a(a)x + a (8a) 

y = c(a)x (8b) 


where a and c are scalar polynomial functions of the parameter a. When a is a specified function 
of time, the system equations constitute a linear time-varying system. We have ( for zero initial 
condition): 

y(t) = c(a(t)) [ e X a ( a ( s )) ds a(T)dT 

Jo 


whose kernel is: 


e 


f)a(a(s))d 


s 


which cannot be written as in Example II.l. Since the state space matrices are functions of a, we 
can describe a- dependent phenomenon with this model. In Section IV, the model given by equations 
(8) is used for aircraft parameter estimation. 
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III. Algorithms for Model Parameter Estimation 


A number of techniques that go under the name of subspace algorithms i6, ^ have been developed for 
unstructured LTI model parameter estimation. In this section, we describe residue-based algorithms 
for model parameter estimation of structured LTI and LPV models. These algorithms have finite 
data guarantees for noise-free data and asymptotic guarantees for noisy data even in the multi-input 
multi-output (MIMO) case provided that the system that generated the data is in the model set. 

Let {uk}k =1 denote the applied input and \yk}k-i denote the measured output. Fix a model 
structure. A model is then defined by a model parameter vector that belongs to some set V. 
Choose a model by fixing model parameter values to p € V. Let us denote by: 



the output obtained by driving the chosen model with the input {rq, ■ ■ ■ , un}. The parameter 
estimation problem is to find a model parameter that minimizes: 

T7 5Z (Vk ~ VkY (Vk ~ vl ) 

k = l 

This is not a standard least squares problem because some of the parameters (e.g. model order) may 
be discrete-valued and, more generally, constrained. To reduce the computational requirements, 
the model order (and all other discrete- valued parameters) are fixed and the best fixed-order model 
is calculated. The model order and other discrete-valued parameters are then selected using infor- 
mation criteria such as those of Akaike (AIC), Bayes (BIC) and Rissanen (MDL) by performing 
a sequence of fixed order model identifications. Our focus in this section is on fixed order model 
identification. 

Consider the fixed order structured LTI model given by equations (1) and (2). The model param- 
eters are the parameters in the state space matrices A, B. C and D. The parameter estimation 
problem can be written as: 

1 N i « 

min — 'V' iyk — y k )' (y^ — y p k ) = min min — ( y & — y k Y (y^ — y k ) (9) 

P eV N ^ [ y K (R*,»oCd)eA b,c,d N ^ Kyk yk> ykl 

where A specifies the apriori information on the parameters in A. We will now show that the 
inner- most minimization problem over the matrices B. C and D is a linear least squares problem. 

Suppose that the pole locations are fixed as in the inner minimization problem. Then, due to its 
modal structure, the A- matrix is fixed by its definition in (2). Choose any B , C and D. Then, we 
have: 


y(s) = G(s)u(s) = [d + C (si — A) 1 B^j u(s) 
= (d + RC (a/ - a) _1 Bj u{s) 

6 of 18 


American Institute of Aeronautics and Astronautics 


u{s) 


= [ D 
= [ D 


R 


C (si - 1) 
R] T(s)u(s ) 


-l 


B 


( 10 ) 


where A, B and C are known matrices independent of B. C and D. R is the block-row matrix 
of residues that depends on B and (7, and the definition of the intermediate transfer function T 
is clear from (10). Note that T is a known transfer function. Figure 2 shows the relationship in 
equation (10) using block diagrams. 


u 


G 


y 



Figure 2. An equivalent way of generating output from input 


The main observation is that the intermediate signal w defined as: 

w(s) = T(s)u(s) 

(see Figure 2) is known since both the applied control input u and the intermediate transfer function 
T are known. Thus, the inner minimization problem can be restated as: 

1 N 1 N 

bcdN^ (Vk - VkY (vk - Vk) = R ] w k)'{yk-[D R]wk) (11) 

’ ’ k=l ’ k=l 

which shows that it is a standard linear least squares problem. A step by step procedure is given 
below. 


Algorithm III.l (Residue Algorithm) Let ({j"*}, {a;,, £,}) 6 A be the given pole locations and 
A be the associated modal matrix as in (2). 

Step 1 : Form the matrices A, B and C . 

Step 2 : Obtain the intermediate signal w by calculating the response of the intermediate transfer 
function T to the input u. This can be done by forming a state space model for T using A, B and 
C , and integrating the state equations. 

Step 3 : Solve the least squares problem (11) to get D and R. 

Step 4: Reduce (A, B , RC. D) to get a state space model that best approximates the data. 


Thus, an algorithm to determine the unknown parameters of a structured LTI model is the following: 
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Algorithm III. 2 (Structured Identification Algorithm) The following steps identify a struc- 
tured LT I state space model given by equations (1) and (2) from input-output data. 

Step 1 : Initialize the parameters of A matrix ( pole locations ({r^}, {oij, Ci})J 

Step 2: Loop until convergence 

(a) Apply residue algorithm to solve the inner minimization problem. Determine the least squares 
approximation error from Step 3 of the residue algorithm. 

(b) Calculate (numerical) derivative of least squares approximation error with respect to pole 
locations. 

(c) Update pole locations using Newton’s method and return to step (a). 

Step 3: Upon convergence, form state space model using the residue algorithm. 


For the LPV model given by (7), a similar algorithm can be given provided that the A-matrix is 
diagonal: 

A(a) = diagonal {r 1 (a), ■ ■ ■ ,r n (a)} (12) 

In this case, algebraic manipulations similar to those used in the derivation of residue algorithm 
show that the LPV model (7) is equivalent to a model of the form: 

x = A(a)x + Bu (13a) 

y = C(a)x + D(a)u (13b) 

where B is a known constant matrix. As before, by moving the unknown parameters in the B- 
matrix of the original LPV model (7) into C matrix of the modified LPV model (13), the overall 
identification problem can be cast as an outer minimization problem over the parameters in A and 
an inner least squares minimization problem over the parameters in C and D. Thus, we have the 
following procedure for LPV model identification. 

Algorithm III. 3 (LPV Model Identification Algorithm) The following steps identify a LPV 
state space model given by equations (7) and (12) from input-output data. 

Step 1 : Initialize parameters of A matrix 

Step 2: Loop until convergence 

(a) Construct the modified LPV model (13) and simulate equation (13a) to get state time histories. 

(b) Formulate and solve least square minimization problem to get the parameters in C and D. 
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(c) Determine the least squares approximation error. 

(d) Calculate (numerical) derivative of least squares approximation error with respect to param- 
eters of A matrix. 

(c) Update parameters of A matrix using Newton’s method and return to step (a). 

Step 3: Upon convergence, form LPV state space model from steps 2a- 2b. 


IV. F-16XL Data Analysis 


This section presents an overview of the F-16XL data and the identification of unsteady aerody- 
namic models from the data using the methods described earlier. 


A. F-16XL Wind-Tunnel Data 

The F-16XL wind tunnel test data used in this study was collected from static and dynamic tests 
performed as part of other work ’ 10,12 using the NASA Langley 12-Foot Low-Speed Wind Tunnel. 
A sketch of the 10%-scale F-16XL model is shown in Figure 3. For the tests, the model was mounted 
on a dynamic hydraulically actuated test rig through a six-component strain gauge balance. The 
mounting was capable of rotating the model through an angle of attack range of —5° to 80° about 
the center of gravity location at 0.558c. The maximum capability of the dynamic test rig was 
260 deg/sec in pitch rate and 2290 deg/sec 2 in pitch acceleration. The tests were conducted at a 
dynamic pressure of 4 psf and a Reynolds number of 10 6 based on the mean aerodynamic chord. 
The data were collected with the leading-edge flaps at 0° deflection. Data were sampled at 100 Hz 
and subsequently filtered using a 6 Hz low pass filter. 

The test data set consists of time histories of the lift coefficient Cl and the pitching moment 
coefficient Cm from static experiments and dynamic experiments with small amplitude oscillatory, 
large amplitude oscillatory and ramp test inputs. Static data were obtained for angles of attack 
from —5° to 80° at zero side-slip and zero trailing-edge surface (flaperons and elevons) deflection. 
The angle of attack (a) test input used to collect the oscillatory data has the form: 

aft) = ao + a a sin(27r/f) 

where ao is the mean angle of attack, a a is the amplitude and / is the frequency in Hz. 

For small amplitude experiments, 

a A = 5°, 20.8° < a 0 < 61.1° and / = 0.6, 1.0, 1.41, 1.75 and 2.94 

Thus, each data file corresponds to an experimental condition indexed by (ao,/) and contains 
information related to system behavior in the neighborhood of (ao,/). The large amplitude data 
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(a) F-16XL 10% scale model dimensions 


(b) F-16XL 10% model in wind tunnel 


Figure 3. F-16XL scaled model used for data collection 


set contains 18 data files of Cl and Cm with 

a 0 = 35°, a a = 10°, 20° and 35° and / = 0.083, 0.25, 0.42, 0.58, 0.75 and 1 

The ramp data consist of several input-output pairs, the input being angle of attack a and the 
output being Cl and Cm- Each pair is from a different experiment in which a was varied as: 

a(t) = otQ + at 


where 


Record 

a (deg/sec) 

a 0 (deg) 

Record 1-9 (ramp up) 

18, 55, 92, 128, 165, 220, 1, 5, 10 

0 

Record 10-15 (ramp down) 

-18,-55,-92,-128,-165, -220 

70 


Figure 4 presents plots of CL{a(t)) from small amplitude oscillatory test data set at two frequencies 
and three values of ao- At ao = 20°, the phase difference between a and Cl is essentially constant 
for both frequencies. At higher ao, however, the variation of Cl with a changes with ao and 
frequency. This behavior can be explained by an indicial response function model whose parameters 
depend on a and a. Then, the model can be reformulated to the state-variable form given by 
equations (6). 
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(a) / = 0.6 rad/s (b) / = 2.94 rad/s 

Figure 4. Small amplitude oscillatory Cl data at different testing conditions and inputs 




(a) / = 0.25 Hz, a A = 35° (b) / = 1 Hz, a A = 35° 

Figure 5. Large amplitude oscillatory Cl data at different testing conditions and inputs 
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B. Model Structure Development 


For one degree of freedom motion in pitch, the fundamental relations for drag, lift and pitching 
moment have the form: 

C' a (t) = C a ( a(t),q(t )) 

where a = D, L or m. Tobak and Schiff formulated these aerodynamic coefficients as 

C a {t) = C a ( 0) + J C aa ( t - t; o(t), 9(t)) a{r)dT + y C aq (t- r; a(r), g(r)) g(T)dT (14) 

where C a (0) is the initial condition and C (lo and C a<j are indicial functions representing the change 
in aerodynamic coefficient due to (separately applied) unit steps in a and q. Following Klein and 
Noderer 8 and Klein and Murphy, 1 we write these indicial functions as: 

Ca a (* - t; a{r),q{r)) = C a<x (oo; a(r), q(r)) - F aoi ( t - r; a(r), q{r)) 

C aq (t ~ T-,a(T),q(r)) = C Uq (oo; a(r),q(r)) - F a<] ( t - T;a(r),g(r)) 

where the deficiency functions F (lo and F aq tend to zero as t — r tend to infinity. Substituting these 
equations into (14) and assuming that the initial condition C a (0) is zero, we obtain: 

C a (t ) = C a (oo; a(t),q(t)) - J F Uo (t - r; a(r),g(r)) a(r)dT ~ y j Q F a q (i - r; a(r),g(r)) q(r)dT 

where C a is the total (due to a and q) aerodynamic coefficient that corresponds to steady flow with 
both a and q fixed at their respective instantaneous values. This is the indicial response function 
model of aerodynamic coefficients given by Klein and Murphy. 10 Reference 10 goes on to suggest 
expressing the deficiency functions around q = 0 as: 

N 

F(t;a, 0) = J2 c i( a ) e ~ ait 

i — 1 

and expanding C a in Taylor’s series around q = 0. This leads us to the model structure for lift 
coefficient 

C L (t) = C L (a(t)) + rj-C L (a(t)) a(t) [ c ai (a(r)) e _a “ i(t_T) d(r)dT (15) 

1 M r* 

- t 7 / c <z* ( a ( T )) e- aqi{t ~ F a(T)dT 

v f~i Jo 

where we used the fact that, in one degree of freedom pitch motion, q = a. The model structures 
for drag and pitching moment are similar. 

Define the state space matrices 

A — diagonal {o a i , ■ ■ ■ , o, a jy , , ■ ■ ■ , dqM } j F — 
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C*(o) [pal (tl) 7 " ' ' 5 (^0 "> Cq\ (ct) 5 ' ' ' 5 CqM (tl)] 

and note that (15) is equivalent to the structured state space model 


i(t) = Ax(t) + B [“g 

y(t) = C(a(t))x 

C L (t) = C L ( a{t )) + yC Lq ( a{t )) a{t) + y{t) 

where the A-matrix is independent of a. More generally, we may take the diagonal entries of A to 
be functions of a to create an LPV system. With N = 1 and M = 1, such a model takes the form: 


x{t) = 


x(t) + 


a{t) 

.«(*). 


' a\ ( a(t )) 0 

0 a 2 {a(t)) 

y{t) = [a {a{t)) c 2 {a{t))]x 

C L (t) = C L ( a(t )) + yC Lq ( a(t )) a{t) + y(t) 


(16a) 

(16b) 

(16c) 


which will be used for identification. The a-dependent terms ai, a- 2 , c\, c 2 , CL(a(t)) and CL q (a(t)) 
are all assumed to be polynomials in a. 


C. Identification Results 

Though it is possible to estimate all the model parameters simultaneously, we follow the sequential 
approach given below to estimate model parameters: 

1. Determine polynomial coefficients of the static term C^(a(t)) using static test data. 

2. Adjust the small amplitude oscillatory aerodynamic coefficient data by subtracting out the 
contribution from static term due to oscillatory a input: 

C L (t) = C L (t) - C L (a(t )) 

where Ci(t) denotes the measured lift coefficient in response to small amplitude a input, and 
CL(a(t)) is the identified static term at the test input a(t). 

3. Determine polynomial coefficients of the damping term from the adjusted small amplitude 
test data 

4. Adjust the large amplitude oscillatory aerodynamic coefficient data by subtracting out the 
contribution from static and damping terms due to oscillatory a input: 

C L {t) = C L (t) - C L (a(t )) - yC Lq ( a(t )) d(i) 

where Ci(t) denotes the measured lift coefficient in response to large amplitude a input. 
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5. Identify a LPV state space model of the form (16a- 16b) to best fit the adjusted large amplitude 
oscillatory data using the residue-based algorithm for LPV models. 


We obtained: 

C L (a) = -0.004 + 1.72a + 11.8a 2 - 32.2a 3 + 26.77a 4 - 7.4a 5 (17) 

from static data. The small amplitude oscillatory data corresponding to the frequency / = 1 Hz 
were adjusted for identifying Ci q . This resulted in 


i-C Lq = -0.7 + 2.86a - 2.05a 2 


(18) 


Finally, the large amplitude data corresponding to cxa = 35° and / = 0.75, 1 Hz were adjusted and 
used to identify the following LPV model: 


x{t) 

v{t) 


-9.7 0 

0 - 8.8 


[ci(a(f)) c 2 


x(t) + 

T 

.0 

O' 

1 

'«(*)' 

.«(*). 

(«(*))]* 





(19a) 

(19b) 


where 


Ci(a) = 0.88 + 6.55a - 5.12a 2 and 02 (a) = -2.03 - 1.07a + 8.88a 2 - 5.53a 3 (20) 


The dependency of A matrix on a was found to be small and neglected in (19a). The complete 
model for Cl is given by (16c). Figures 6-7 show comparison of response of identified model with 
the data used for identification and with the large amplitude oscillatory data at a a = 20° . 




(a) ola = 35°,/ = 0.75 Hz 


(b) ola = 35°, / = 1Hz 


Figure 6. Comparison of response of identified model and large amplitude test data 


The ramp data was used to test the validity of the identified model. Figures 8-9 show comparisons 
of experimental data and output of identified model. These figures show good agreement up to a 
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(a) q a = 20°,/ = 0.7572* 


(b) a A = 20°,/ = 1 Hz 


Figure 7. Comparison of predicted response of identified model and large amplitude test data 


ramp rate of 18° / s. A part of the error is due to initial conditions which where set to zero for the 
identified model. The largest percentage error occurred at a ramp rate of 18 ° /s as shown on the 
right in Figure 9. 


V. Conclusions and Directions for Future Research 


This paper presented local and parameter dependent model structures for modeling unsteady aero- 
dynamic effects. We have also developed additional nonlinear model structures details of which 
can be obtained from the authors. Table 1 gives a summary of the model structures. The unstruc- 
tured model identification can be carried out using the subspace method and its extensions to the 
nonlinear case. The structured model identification can be carried out using the residue algorithm 
and its extensions. 

The algorithms were evaluated on a limited set of wind tunnel data for the F-16XL provided 
by NASA Langley Research Center. Our experience has been that the model structures and 
algorithms are well suited for modeling unsteady and nonlinear aerodynamic effects, in particular 
the structured models. It should be noted that the key issue of model structure determination 
is still un-resolved. We have no systematic procedure to determine model structures (other than 
trying out all the structures in a library of model structures). This issue is likely to remain difficult. 
The identification results presented in this report dealt with an aircraft with one lifting surface. 
This approach should be extended to a wing-tail configuration (some work along these lines can be 
found in Reference 11). In addition, modeling should be extended to aircraft non-planar motion 
taking into account the effect of angle-of-attack and side-slip on parameters. One of the main issues, 
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(a) Ramp up from oto = 0° at rate a = 1° / s (b) Ramp up from qo = 0° at rate a = 5°/s 

Figure 8. Comparison of predicted response of identified model and validation ramp test data 




(a) Ramp up from qo = 0° at rate a = 10° / s (b) Ramp up from ao = 0° at rate a = 18° / s 

Figure 9. Comparison of predicted response of identified model and validation ramp test data 
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Property 

Unstructured 

LTI 

Structured 

LTI 

Hammerstein 
and Wiener 

LPV 

Parameter-Dependent 
Hammerstein /Wiener 

Linear 

Yes 

Yes 

No 

Yes 

No 

Global 

No 

No 

No 

Yes 

Yes 

Apriori 

Information 

No 

Yes 

Yes 

Yes 

Yes 

Algorithm 

Subspace 

SID 

Modified SID 

Modified SID 

Modified SID 


Table 1. Summary of Model Structures (SID-Structured Identification Algorithm) 


as with the problems studied in this report, lies in model structure determination and analytical 
studies to gain a physical understanding of the flow effects should be undertaken. We suggest three 
additional areas for future research. 

Identification of Reduced-order Models of CFD Solvers: Unsteady CFD solvers make the 
fewest assumptions about the flow and, in many cases, they can accurately predict flow around 
geometries. But, their inclusion in aircraft analysis and design is not feasible due to the high com- 
putational cost. The development of computationally efficient reduced order models of unsteady 
flow is thus of paramount importance. Though the focus of our work has been modeling of aerody- 
namic effects for flight dynamics, the models and algorithms described in this paper can be applied 
to obtain reduced order models of CFD solvers. It has potential applications in aeroservoelasticity 
and control of morphing wings. 

Uncertainty Description: Nonlinear phenomenon such as Hopf bifurcations are properties of 
mathematical models that occur at specific parameter values. If these properties persist with 
respect to uncertainties, then we can be reasonably confident that the true system will also exhibit 
such properties. So, we suggest that useful models be required to supply a description of the 
uncertainties so that phenomena observed during analysis of an identified model can classified as 
artifacts of the model or physical phenomena. Computation of modeling uncertainties from data is a 
difficult problem. The maximum likelihood estimation (MLE) and Bayesian parameter estimation 
techniques may be able to determine parametric uncertainties. A deeper look into this issue is 
suggested for future work. 

Experiment Design: The reliability of any model identified from experimental data is dependent 
to a great extent on the experiment itself. The problem of input design, also known as experiment 
design, has been the subject of vigorous research. At present, experiment design for linear system 
identification is well understood and there is a general framework for addressing constrained and 
nonlinear systems. See for instance the references 7, 14, 18, 19. However, experiment design for 
high angle of attack flight dynamics and unstable flight regimes has not received the same level of 
attention. The use of experiment design methods in industry can reduce the amount of testing and 
improve safety. 
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